{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# GP Regression with Uncertain Inputs\n",
    "\n",
    "## Introduction\n",
    "\n",
    "In this notebook, we're going to demonstrate one way of dealing with uncertainty in our training data. Let's say that we're collecting training data that models the following function.\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "y &= \\sin(2\\pi x) + \\epsilon \\\\\n",
    "  \\epsilon &\\sim \\mathcal{N}(0, 0.2) \n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "However, now assume that we're a bit uncertain about our features. In particular, we're going to assume that every `x_i` value is not a point but a distribution instead. E.g.\n",
    "\n",
    "$$ x_i \\sim \\mathcal{N}(\\mu_i, \\sigma_i). $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Using stochastic variational inference to deal with uncertain inputs\n",
    "\n",
    "To deal with this uncertainty, we'll use variational inference (VI) in conjunction with stochastic optimization. At every optimization iteration, we'll draw a sample `x_i` from the input distribution. The objective function (ELBO) that we compute will be an unbiased estimate of the true ELBO, and so a stochastic optimizer like Adam should converge to the true ELBO (or at least a local minimum of it)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The autoreload extension is already loaded. To reload it, use:\n",
      "  %reload_ext autoreload\n"
     ]
    }
   ],
   "source": [
    "import math\n",
    "import torch\n",
    "import tqdm\n",
    "import gpytorch\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "%matplotlib inline\n",
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Set up training data\n",
    "\n",
    "In the next cell, we set up the training data for this example. We'll be using 20 regularly spaced points on [0,1].\n",
    "We'll represent each of the training points $x_i$ by their mean $\\mu_i$ and variance $\\sigma_i$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Training data is 100 points in [0,1] inclusive regularly spaced\n",
    "train_x_mean = torch.linspace(0, 1, 20)\n",
    "# We'll assume the variance shrinks the closer we get to 1\n",
    "train_x_stdv = torch.linspace(0.03, 0.01, 20)\n",
    "\n",
    "# True function is sin(2*pi*x) with Gaussian noise\n",
    "train_y = torch.sin(train_x_mean * (2 * math.pi)) + torch.randn(train_x_mean.size()) * 0.2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x12099f470>"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeQAAADBCAYAAAAeuMPIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAUoElEQVR4nO3df2zc9X3H8ZeviXMYsjQJ9VwFutJAmEKbSeiNECRUrnBXFYRUplaFppHG4nNZhbUWupCSrFamNA4sa2cDXXutVZqk60y1AqqKhjRFzAWPlncEFZRUSRmN2Tj7UJIrVOTiZPb+8DlcHPt8vl/fz9nPh3TS977fu2/efDjf6z7fH59Pw/j4uAAAQLRiURcAAAAIZAAAgkAgAwAQAAIZAIAAEMgAAASAQAYAIACLovzHt27dyj1XAIAFZ/fu3Q1T10UayJK0Y8eOqEsoSjqdVnNzc9Rl1AXaqni01dzQXsWjrYpX67bq6uqadj2HrAEACACBDABAACI/ZI35LZVKadOmTdq/f79aWlqiLgdADZ05c0bDw8M6deqUQh6meWxsTG+99VbF99vQ0KAlS5aopaVFixbNHrcEMqqqu7tbg4OD2rVrl3p7e8veHwEP1I/h4WFdeOGFWrVqlRoazruGKRinT5/W4sWLK77f8fFxZTIZDQ8P65JLLpn19QQyzjEZeM8880xF95tMJpVMJiuyr1gsVrGAB1A9p06dCj6Mq6mhoUHvfe97dfz48aJezzlknGOyRxuysbExJZNJxePxsh8Aqmd8fLzoME6lUmpra9Pw8HCVq6qthoaGog/X00OeBwgWAPWu3NNb27dv15VXXqk9e/boK1/5il566SU98MADM76+q6trxttuR0dHdd999+nEiRO68cYb9dxzz+mWW27Rxz/+8TnXNRcEMmriC1/4goaHh9Xf31/yPjo7O9XX16fGxkaNjo6qvb2dw9ZAHZuuM5F/eiubzRa9r7vvvlsrVqzQ3r17tWnTJp04caLg6wuNgdHY2Kh169ZJkj73uc/p05/+tDZs2KCrr75aK1euPPu6bDarffv2KZFIFF1nISUHspl9VFKXu984Zf31ktZr4nD49909XV6JmM1cPrSFVCrwqnWTfTqdViKR0ObNm9XX1zfvDm0BKN2KFSvOeb58+XLt2LFD4+Pjampq0ttvv621a9fq97//vW6++WZt375du3fv1he/+EV94hOf0IEDB/Too49Ou+/GxkatX79eP//5z3Xo0CFdddVVevHFF3XbbbfpJz/5iVpbW7V3796z+7/zzjtL+m8o+Ryyuw9IumCaTd2S9kj6kaT6GIYLkt4NvIGBASUSCY2MjERd0jn6+/vV09OjdevWqaenp6zeNoDoZbNZZbNZJRIJxWIxxeNxxWIxdXR0VKSjcdlll+mGG27Qli1b9KlPfUrNzc06ePCgLr30Ui1evFjvf//7tWLFCm3evFmZTGbW/cXjcd18881atmyZjhw5ojVr1mjVqlW64oorztl/qcq9qGs0/4mZrZF0xt3H3X1I0g1l7h81ROABiEI1OwOLFi3S+Pi49u7dq9bW1vMusJq86KzQhVenT5/WCy+8oOuuu04/+MEPdM0112jJkiVntxfa/5xqLfmd07tYUv6B++WzvSGdro8j2sX8esIE2qp4tNXc0F7FC6GtxsbGdPr06Vlft3///rPLe/bskaSi3jfVkSNH9Nprr+k3v/mNVq9ercOHD+v48eO67rrr9Otf/1qPPPKI0um03njjDR05ckSvvvqqXn31Vf3ud7/TsWPHzr5vdHRUBw8e1MmTJ3XRRRfpV7/6le6//341NTVpaGhIfX19SqfTOnr0qDKZjF566aVz9p9Op7V8+bvxNzY2VlTWNZST5mb2tLu35j2/QlKPu9+Ue/7f7v6hmd6/devWcSaXmH9oq+LRVnNDexUvhLb67W9/q8svvzzSGopRrYFBJk1th66ururN9mRm75HU5O5HzCyeW/chSU9XYv8AAMx35Vxl/RFJq83sw5Iuk7RB0r2SdprZFklxSdsqUiUAAPNcyYHs7i9JujT39GVJP82tPyDpQPmlAQDq2eQoVQt16ExpbqOVMXQmAKAqlixZokwmE/RMT9U0OblE/hXZhTBSFwCgKlpaWjQ8PKzjx48HHcpjY2OKxSrfP82ffrEYBDIAoCoWLVpU1LSDUQvhinSJQ9YAAASBQJ4H5uu0ZQCwkBDI80D+tGUAgPrEOeQIVGv+4vxpyyqhUrNIAQBmRw8ZAIAA0EOOQCV7npWawxgAEC16yHUu9DmMAQDFoYdc5/LnLO7p6YmwEgBAOeghAwAQAAIZAIAAEMhAoBjwBVhYCGQgUAz4AiwsXNQFVFElBoEpd8AXBngB6gM9ZAAAAkAgA1WUzWYLPoaGhqZdn0gkFIvFFI/HFYvF1NHRMeu+ZnoAqA8EMhAgBnwBFh7OIQMBYsAXYOGhhwwAQAAIZAAAAkAgAwAQgJID2czuMbNNZnbXNNseN7NhM/tueeUB1cVoWABCUVIgm9kGSSvdfZ+k5WZ2bd62ayT9s7u3uHuiQnUCVcFoWABCUepV1jdJOpRbfiX3/Be55x+T1GlmByT9tbu/U16JwPQqMQrWpHJHw5IYEQtAeUoN5IslncgtZyW1TG5w9wfM7BuS7pe0VdLXCu0onU6XWEJtZTKZqEuoGwu1rUr5LC/UtioV7VU82qp4obRVqYH8pqSm3PJSScfyN7r7GTO7V9L3Z9tRc3NziSXUXj3VGrVatFW5PdLOzk719fWpsbFRo6Ojam9vV29vb4WqKx6fq7mhvYpHWxUvhLYq9aKuJyWtyy2vlfSUmS2TJDNryK1fKumZ8soDqofRsACEpKQesrs/a2YfM7M7JGVyj29Lul3SM2b2gqQXJH2vYpUCFcZoWABCUvLQme6+c8qq23Pr15dVEQAACxADgwAAEAACOQIMRgEAmIpAjgCDUQAApmL6xVlUcvCJqSoxGMVUDE4BAPWJHjJQQSGfjgi5NgAE8qyy2ayy2ayGhobOLpfzSCQSisViisfjisVi6ujoqMh+Jx+IVsinI0KuDQCHrGtucjCKzZs3q6+vj95KYCp1iqLc0xGFflyVW2M5tfGjD6geArnGGIwCADAdDlkDeWp9OmKmUyGVrrFSp0oAVA+BDFRIyGNjh1wbgAkcsgYqJOTTESHXBmACPWQAAAJAIAMAEAACGQCAABDIAAAEgEAGACAABDIAAAEgkAEACACBDABAAAhkAAACQCAXIZVK6TOf+QwzMwEAqoZALkJ3d7d++ctfMo8sAKBq5t1Y1pWaz3Y65c5xOxNm0QEAlBzIZnaPpLSkZe7+UN76NZI+K+kdST9198NlVwkAQBWkUinddttt6u/vV0tLS6S1lHTI2sw2SFrp7vskLTeza/M290j6pqSHJO0uv8S5KWc+20LzyC5ZsqSseWSZYxYAwhPSKclSe8g3STqUW34l9/wXZnaBpNXu/gdJMrPLzGyRu58pv9RoTM4je+utt+qxxx7jwi4AqDPFnMqc6ZRkLTtNpQbyxZJO5Jazkib7+cslvZX3ujOS3icpNdOO0ul0iSXUxoMPPihJymQy2rZtm6Twa45aJpOJuoS6QVvNDe1VPNqqMmr5fV9qIL8pqSm3vFTSsdzyMUn5P0WaJBX8VDQ3N5dYQu3VU61Ro62KR1vNDe1VPNpqwky93M7OTvX19Wnx4sU6ffq02tvb1dvbW+Pq3lXqbU9PSlqXW14r6SkzW+bupyQdNbMmM4tLet3dT1aiUAAAKmnylOQTTzyhRCKhkZGRSOspqYfs7s+a2cfM7A5N9IAzkr4t6XZJ90raIumUpLsrVSgAAJXU398vaSKYW1tboy1GZdz25O47p6y6Pbf+ZUkvl1MUAAALDSN1AQDqRiqVUltb27y844VABgDUje7ubg0ODgZx33ClzbuhMwEA9a2c+4al+h2OmB4yAAABIJABAEGZbSjjeDxecCjjekUgAwDqwuR9wwMDA0HcN1xpnEMGANSFyfuGJamnpyfCSqqDHjIAAAEgkAEACACBDABAAAhkAAACQCADABAAAhkAgAAQyAAABIBABgAgAAQyAAABIJABAAgAgQwAQAAIZAAAAkAgAwAQAAIZAIAAEMgASpJKpdTW1qbh4eGoSwHmBQIZQEm6u7s1ODioXbt2RV0KMC8smusbzKxZ0l2ShiW96O6DU7Z/UNKgJsL+8+7+HxWoE0CVxOPxst6fTCaVTCbn/L5sNnvO81QqpU2bNmn//v1qaWkpqyagHpXSQ94lab+7f0vSV82sYcr2z0r6E3dvIYwBFIseNxa6UgL5zyUdyXv+wckFM2vMbT9qZhvLKw1ALWSz2Tk/EomEYrGY4vG4YrGYOjo6ztlejHg8fs4jmUxqbGxMyWTyvG2Tjw984APTrgfmg4KHrM3sPklrpqx+n7uP55azklokvSZJ7j4q6UYzu0TSz8zseXc/XOjfSKfTJRVea5lMJuoS6gZtVbx6bavXX39dGzdu1MaNG/XDH/5QQ0NDkf4t18v3SC3V62crCqG0VcFAdvfzjh2Z2fq8p0slHZvmff9jZl+X9GFJBQO5ubm5uEoDUE+1Ro22Kl49ttXjjz9+drm1tfW87cX2kid1dnaqr69PjY2NGh0dVXt7u3p7e897XTqdrsv2igptVbwQ2qqUQ9ZPm9nlueUl7n7YzJaa2XskKe+c8gWSnqtEkQDmt3Q6rUQioYGBASUSCY2MjERdElBzc77KWlKXpE4zG84tS9LXNBHUxyU9bGY/lvSsu79RoToBzGP9/f1nl3t6eiKsBIjOnAM5F7JfnbLub/OeXl1uUQAALDQMDAIANcYoZ5gOgQxg3gk98LjnGtMhkAHMO1EF3kz3T5dyz/V8udc69B9HISnloi4AqKlSA6mUYT3nessWCsv/cTTdrWx4F4EMABVSTJgXe891vSj2x9JsP474IcQhawB1YHJIzqGhobKH9ZztUW3cc42Z0EMGMK9MBt7mzZvV19cX3LnL+XbPdaEfMfPtaEC1EcgA5pX5Fnj1LPQfR6EhkAEAVcGPo7nhHDIAAAEgkAEACACBDABAAAhkAAACQCADABAAAhkAgAAQyAAABIBABgAgAAQyAAABIJABAAgAgQwAQAAIZACoc6lUSm1tbUzeUOcIZACoc93d3RocHNSuXbuiLgVlYLYnAAhUPB6f0+uTyaSSyWTB1xSavxjRKimQzewCSV+WNObuu6fZfo+ktKRl7v5QeSUCADD/lXTI2t1PSnJJ5/18M7MNkla6+z5Jy83s2vJKBICFKZvNzvpIJBKKxWKKx+OKxWLq6OhQNpvV0NDQtK9HuMo5hzw6w/qbJB3KLb+Sew4AqIJ0Oq1EIqGBgQElEgmNjIxEXRJKNOshazO7T9KaKasfl5SZ4S0XSzqRW85Kaim0/3Q6PVsJQchkZvrPxVS0VfFoq7mhvc734IMPnl3etm2bpInv1XLaamRkRHfddZcefvhhNTc3l11j6EL5XM0ayO4+7WV7ZtY6w1velNSUW14q6Vih/dfT/+x6qjVqtFXxaKu5ob2KV2pb7dy5U88//7ySyaR6e3srXFWYQvhcVewqazNbJuktSU9K+qSkRyWtlfTvlfo3AADlmcuV27Ndtc056coq6RyymS2SdL2kq8xseW71dyR9xN2flZQ1szskZdx9oDKlAgAwf5XUQ3b3M5J2TVl3W97yzjLrAgBUwWy92s7OTvX19amxsVGjo6Nqb29fMIeto8ZIXQCAs7hqOzqM1AUAOKu/v//sck9PT4SVLDz0kAEACACBDABAAAhkAAACQCADABAAAhkAgAAQyAAABIBABgAgAAQyAAABIJABYBapVEptbW0aHh6OuhTMYwQyAMyiu7tbg4OD2rVr2tlogYpg6EwAC9JcpiGcNNt0hFMxPSHmgh4yAAABIJABLEjZbLaoRyKRUCwWUzweVywWU0dHR9HvBeaCQAaAApiOELXCOWQAKIDpCFEr9JABAAgAgQwAQAAIZAAAAkAgAwAQgMgv6urq6oq6BAAAItcwPj4edQ0AACx4HLIGACAABDIAAAEgkAEACACBDABAAAhkAAACEPltTyEys3skpSUtc/eH8tavkfRZSe9I+qm7H46oxGAUaKvbJX1J0h9J2uTuHlGJQZmpvfK290na5+5P17q20BRqKzP7U0k3SHrZ3f8rivpCUuDv8FZJK3NP33H3f4mivpCY2Ucldbn7jVPWXy9pvSY6qt9393Sta6OHPIWZbZC00t33SVpuZtfmbe6R9E1JD0naHUV9IZmprcysQRN//NdK2iNpR4RlBmOWz5bM7BZJF0VSXGAKtZWZXSkp4e7fJYxn/Vz9jbt/z92/J+mvoqkwLO4+IOmCaTZ1a+L76keK6DuLQD7fTZIO5ZZfyT2XmV0gabW7/8HdT0m6zMwW+hGGadvK3cfd/Ync+uclpSKoLUTTtpckmdllmjhidWia9y1EM7aVpF5JR82sJxdGC12htjpoZn9vZibpWzWvLFyj+U9yRz/P5L67hjRx9KXmCOTzXSzpRG45K6klt7xc0lt5rzsj6X01rCtEM7VVvjZJ36hZRWGbtr1yP+w+6e6PRVVYgGZqqwslfVATR6n+UdKPzawxigIDUujv8O8krZb0D5IGalxXPclvQ2ni+77mCOTzvSmpKbe8VNKx3PIxSfG81zVJytSwrhDN1FaSJDO7XNJRd3+l1oUFaqb2+qikz5vZ05L+UtI/mdmqmlcXlpnaqlHSSXcfy/Vk3tD0PwQXkkJ/h1+XdKcmDsX+a43rqif5bShJp6IogkA+35OS1uWW10p6ysyW5Q5THzWzJjOLS3rd3U9GVmUYpm0rSTKzP5b0Z+7+b2Z2Ua5ns9DN9Nk64O7Xu3urpEckfcnd/zeiGkMxU1udkHTKzCbPtb8pibaa4e9Q0jp3f9vdfyZpcSTVBczM3mNmS939iHIdLjP7kKSno6iHQJ7C3Z+VlDWzOzTRA85I+nZu872Stkj6sqS7o6kwHDO1lZmtlPSUpK+amUv6T01cmb6gzfLZQp5Z2uouSV25K/nvd/f/i6jMIMzSVt80s04z+wtJ34mqxpCY2UckrTazD2vifPv23KadZrZF0uclbYuiNiaXAAAgAPSQAQAIAIEMAEAACGQAAAJAIAMAEAACGQCAABDIAAAEgEAGACAABDIAAAH4fyfRcxByOgJbAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 576x216 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "f, ax = plt.subplots(1, 1, figsize=(8, 3))\n",
    "ax.errorbar(train_x_mean, train_y, xerr=(train_x_stdv * 2), fmt=\"k*\", label=\"Train Data\")\n",
    "ax.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setting up the model\n",
    "\n",
    "Since we're performing VI to deal with the feature uncertainty, we'll be using a `~gpytorch.models.ApproximateGP`. Similar to the [SVGP example](./SVGP_Regression_CUDA.ipynb), we'll use a `VariationalStrategy` and a `CholeskyVariationalDistribution` to define our posterior approximation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "from gpytorch.models import ApproximateGP\n",
    "from gpytorch.variational import CholeskyVariationalDistribution\n",
    "from gpytorch.variational import VariationalStrategy\n",
    "\n",
    "\n",
    "class GPModel(ApproximateGP):\n",
    "    def __init__(self, inducing_points):\n",
    "        variational_distribution = CholeskyVariationalDistribution(inducing_points.size(0))\n",
    "        variational_strategy = VariationalStrategy(self, inducing_points, variational_distribution, learn_inducing_locations=True)\n",
    "        super(GPModel, self).__init__(variational_strategy)\n",
    "        self.mean_module = gpytorch.means.ConstantMean()\n",
    "        self.covar_module = gpytorch.kernels.ScaleKernel(gpytorch.kernels.RBFKernel())\n",
    "        \n",
    "    def forward(self, x):\n",
    "        mean_x = self.mean_module(x)\n",
    "        covar_x = self.covar_module(x)\n",
    "        return gpytorch.distributions.MultivariateNormal(mean_x, covar_x)\n",
    "\n",
    "\n",
    "inducing_points = torch.randn(10, 1)\n",
    "model = GPModel(inducing_points=inducing_points)\n",
    "likelihood = gpytorch.likelihoods.GaussianLikelihood()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Training the model with uncertain features\n",
    "\n",
    "The training iteration should look pretty similar to the SVGP example -- where we optimize the variational parameters and model hyperparameters. The key difference is that, at every iteration, we will draw samples from our features distribution (since we don't have point measurements of our features).\n",
    "\n",
    "```python\n",
    "# Inside the training iteration...\n",
    "train_x_sample = torch.distributions.Normal(train_x_mean, train_x_stdv).rsample()\n",
    "# Rest of training iteration...\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "70396890e54b44979d857ae566d116e4",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=400), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "# this is for running the notebook in our testing framework\n",
    "import os\n",
    "smoke_test = ('CI' in os.environ)\n",
    "training_iter = 2 if smoke_test else 400\n",
    "\n",
    "\n",
    "model.train()\n",
    "likelihood.train()\n",
    "\n",
    "optimizer = torch.optim.Adam([\n",
    "    {'params': model.parameters()},\n",
    "    {'params': likelihood.parameters()},\n",
    "], lr=0.01)\n",
    "\n",
    "# Our loss object. We're using the VariationalELBO\n",
    "mll = gpytorch.mlls.VariationalELBO(likelihood, model, num_data=train_y.size(0))\n",
    "\n",
    "iterator = tqdm.notebook.tqdm(range(training_iter))\n",
    "for i in iterator:\n",
    "    # First thing: draw a sample set of features from our distribution\n",
    "    train_x_sample = torch.distributions.Normal(train_x_mean, train_x_stdv).rsample()\n",
    "    \n",
    "    # Now do the rest of the training loop\n",
    "    optimizer.zero_grad()\n",
    "    output = model(train_x_sample)\n",
    "    loss = -mll(output, train_y)\n",
    "    iterator.set_postfix(loss=loss.item())\n",
    "    loss.backward()\n",
    "    optimizer.step()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAd0AAADFCAYAAAARzygsAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3deXxc5Xno8d/su2ZGu2xLXuR9AQPHNgQCbgMYCFlvIM1CQ0JweglpyG1LaNqU0HtLSZM2DSWhMeW2gQRKgAQuWwAnLAGM7WODsWVbXrXYljSjZUazb2fuHyPLuy2NRjMj6fl+PvOZmaOZc1690sxz3ve87/PqstksQgghhBh/+lIXQAghhJgqJOgKIYQQRSJBVwghhCgSCbpCCCFEkUjQFUIIIYpEgq4QQghRJBJ0hRBCiCIx5vMmRVG8wI8ABfjfqqo+UdBSCSGEEJNQvi3dWuArwNXAZwtXHCGEEGLy0o0lI5WiKKsAm6qqrxesREIIIcQklVf3MoCiKHOA+wA/8PqZXnfXXXdJnkkhhBBTyn333ac73fa8g66qqgcURfkI8IGiKDWqqvrP9Np77rkn38OMmc/no7a2tmTHn6ik3vIj9ZYfqbf8SL3lZ7zr7e677z7jz8Y0ellVVQ14F+gfy36EEEKIqSDf0ct3AIuBt4GfqaqaKWiphBBCiEkor6Crquq/FrogQgghSiedTtPd3U0ikWCyL/mqaRqDg4MF2ZdOp6OmpoaKiooRvT7va7pCCCEmj+7ubhwOB9OnT0enO+0YoEkjlUphMpkKsq94PE5PT8+Ig65kpBJCCEEikcDj8Uz6gFtoFouFdDo94tdL0BVCCEE2m5WAmwedTjeq7njpXhZCCFF0Dz30EDU1NUSjUfR6PX/yJ3/C3//937N69Wouv/zygh9vYGCAr33ta/zqV78CIBwOs3btWpqamli6dClvvfUWa9eu5cILLyz4sY8nQVcIIcQwq9VSsH3F44nTbn/iiVy6/k9+8pMA/O3f/i1vv/02s2bNGrdBXF6v94Trrk6nk0WLFnH55ZdzxRVXcO2117JmzRreffddjMZjobGvr4+XXnqJL37xiwUphwRdIYQQRfXcc8/x7W9/e/j5ihUreP7551m0aBGvvvoqP/rRj7jhhhuoq6sjGAyyceNG7r33Xh566CECgQALFixg+/btZLNZ7HY7zz77LI888gi9vb10d3ejaRqZTIbXX3+dBx54gPvvv585c+awf//+M5apqqqKpqYmtm/fzq9//WsWL15MMBikubmZV155hcsvv5yf/OQnNDU10djYyMc//vG8fne5piuEEGJYPJ4o2O1MNE07oUVrMpkwGAwAXHXVVTz++OM89NBDhEIhXn75ZW655RbWr19PKBTiggsu4PDhw8yePZsPf/jD3Hnnndxxxx1s2LABv9/PJz7xCZ588klqamrwer384Q9/wGQycf311zNnzpxz/v5Wq5VPfvKT1NbWsmXLFubNm8e8efNoamriE5/4BB6Phw8++CDv+pWgK4QQoqiuu+463n777eHnW7du5WMf+9jwc6vVyrRp07jkkku46aabuPXWW0mlUuh0Oq655hquv/56gOFu4E984hM8++yzw4G7p6eHyy67jNtuu43u7m76+3NJE882UGxgYID+/n7mz5/PI488wurVq08YXBYMBnnuuee49NJLx9QFLt3LQgghiuoLX/gCP/nJT3jmmWdIp9MsWrSIVatWkclkWL9+Pfv27ePOO+/kiSeeoKamhquuuoo1a9awbt06Nm7cyB133MG+ffsIBoNcccUVmM1mzjvvPFatWgXAjTfeyOrVq7n++uv55je/yec//3keeOABenp6OHToEHV1dYTDYVpaWkgmkwSDQTZt2sTDDz+MXq9n586dPProo/j9fvR6PTt27CAQCLBhwwaam5vZvXs38Xgcq9U66t99TEv7jcRdd92VlQUPJh6pt/xIveVH6i0/hay3ffv2MXfu3ILsq9wVMjkGnFp3d9999xlXGZLuZSGEEKJIJOgKIYQQRSJBVwghhCgSCbpCCCFEkUjQFUIIIYpEpgwJIYQ4wb+9dubMTSP1jT9qPuPPHnvsMf7hH/6BHTt2DM+D/fM//3OMRiP33XcfZrN5zMcvV9LSFUIIUVQf+tCHqKmp4eWXXwYgFArR1tbG8uXLJ3XABWnpCiGEKIG1a9fys5/9jGuuuYbnn3+ej370o6TTaR588MHh/MoGg4GBgQE2bdrEvffey5e//GVuuOEGfv7zn/PCCy+csDDBRCEtXSGEEEU3c+ZMjEYje/fuJR6PY7PZSKVSJ+RXXrx4MQsXLqSjowOPx0NlZSWf//znmT59Ot3d3aX+FfKSV9BVFMWlKMqTiqIcUBTlp4UulBBCiMlv7dq13HLLLcPr52YymRPyK7/wwgtUVlZSX18PHMudrNfrSafTJSv3WOTbNr8YuBnIAu8pirJCVdXNBSuVEEKIkjnbIKhC2LRpEz6fj69//eu8+uqrNDY2sm3bNvR6PR988MFwfuVEIsHDDz9MNBrlzTffZP/+/bS3t9PV1cX+/fuZNWvWuJZzPOQVdFVVffXoY0VRdgATs50vhBCi6G688cbhx//0T/8EwD//8z+f8rrLLrvshOevvfYaAC+88MI4lm58jekqtKIoLqBDVdXOs73O5/ON5TBjEggESnbsiUzqLT9Sb/mRestPIetN0zRSqVTB9lfOMplMQfenadqI49xYh37dBPzduV5U6tVDSn38iUrqLT9Sb/mRestPoeptcHCwoCvvlLtC/q56vX7Ef4e8Ry8rivJJ4BlVVUOKotTlux8hhBBiqsh39PJtwI+A/6coygfARwtaKiGEEGISyncg1U8BmSokhBBCjIIkxxBCCDEqXV1dXHnllWNKUPHQQw/x1FNP8fzzz3PTTTed8/WpVIrHH3+c2267jZtvvnl4ezqd5mtf+xptbW15l6WYJl4OLSGEECX1j//4j7zzzjvce++93H///aN+/5NPPkk6neYzn/kMAC6X65zv2b59O8lkkp/+9MROVqPRyIwZM0ZdhlKRlq4QQohzslqtw7d169ahaRrr1q0b3jYaTz/99HAWKoArrriCdevWsX79eu688076+vr41Kc+xWOPPcaaNWtIp9O89tpr7NixgyNHjnDDDTcA8Pjjj/Piiy/y1ltvAfCb3/yGp556ittvv53du3fzuc99jh/+8Id84xvfAOCRRx7h5z//Od///vdpa2vjv/7rv7j99tvZu3dvgWrp3CToCiGEKKpMJnPCnOD29nYOHjzIlVdeidfrZdOmTafkWVYUhWXLljFt2jTcbjfxeJxXX32V6667jksvvRTItaBramrwer3MmDEDh8PBX/7lX9La2squXbuIx+N86Utf4nOf+xzr1q2joaEBRVE4ePBg0X53CbpCCCHOKR6PE4/HufXWW9Hr9VitVvR6PWvXriUej49qX2vWrGH9+vXDzzs7O+nszOVYqq6uxu12nzPPciwW4/Dhw8CxnMw9PT1cdtll3HbbbSfMw81ms2iaxs6dO4efZzIZPB4PN998M/PmzRtV+cdCgq4QQogR8/l83Hrrrbz55pvceuut9PT0jHofX/nKV+jr6+MHP/gBv/zlL9HpdMyaNYtf/OIX6HQ6lixZckqe5R07dtDS0kIwGGTfvn2k02mWL1/O3/zN39De3s6ePXu48cYbWb16NY888ggdHR10dXXR3t7OkSNHqK2tJRwO84UvfIGDBw9y66238q1vfYtbbrmlqGv46rLZ7Lge4K677srec88943qMs/H5fJLpJg9Sb/mResuP1Ft+Cllv+/btY+7cuQXZV7lLpVIFzUh1ct3dfffd3HfffbrTvVZGL4tRO/5E7Wi3jhBCiHOToDsFZbNZYimNWDJDLHXcLZnJbU9l8PcHsB9Kkc5kSWtZ0hmNlJYlncmS0bJkOS7wojv6AAC9Dox6PWaDDpNRj9mgx2TQYTbqMRlyz+1mw3E3Iw6LAbvJgF4vQVwIMXlJ0J2kwvE0g/E0g/EUoaHHoUSaUCxFKJEmo539skIkksCRHtm/x3AAHrrLZCGjZUikgcTIy6xDh9Wkx2kx4rabcFuNuG0m3DYTHpsJh8UgLWshxolOpyObzcpnbJRGW2cSdCcwTcsSiKUYiKYYiCYZiKboj+TuUxmt1MUbtSzZ4Va3P3xqtDbq9VRYjXgdZqqdZqocZmqcZipsU2dlFCHGi8ViIRAI4PF4JPCOQiKRwGgceSiVoDtBJFIZ/OEk/nACfyh3H4im0MZ5IFw5SWsa/dEk/dEk+/3HtpsNeqqcZqqdFqqdZmpdFqodZumqFmIU6uvr6e7upr+/n/EeYFtqmqah1xdm8o5Op6OmpmbEr5egW4biqQw9gwl8oUQu0IYSDManxuLS+UhmNLqCcbqCx+YKmgx6al0WGtxW6ios1FdYsZsNJSylEOVtoqVTHItSjpaXoFtimpalN5KkezBOz2CCnsFcC/b4gUpi9FIZjcOBGIcDseFtbpuJBreV6R4rjV47Lqv8+wshiku+dYoskcpwJBinezDBkWAcfygxIa+/TkTBWIpgLMXu7hCQC8IzvDYavTZmeGzYpCUshBhnEnTHWTie5vBQ1+eRQJz+SFJasWXiaBBuOTKIDh2VDjNNlTZmVtmZ5rZikGvCQogCk6BbYAPRJEcCcQ4H4hwJxgnJtdgJIUuWvkiCvkiC9zoDmA16GivtzK6yM7PKLteDhRAFIbmXxyCbzdIXTvLB4SC/benh/77dzi82dvL7Vj+tPaGSBdzBfj8//asvMdjvP/eLS7C/8dpnISUzGvv9Ydbv9vF/327nV+phNrUN4A+NYuKxEEKcRILuKGSzWXyhBO93BnlxRzcPv93OY5s7eWNPL3t9YSLJU1fCKIX1jz1IW8sW1j/2YFnub7z2OV6yZOkJxdl4sJ//Vg/xyLsdvLO/j57B0a2sIoQQ0r18FpqWC7KHh67HdgXjJNKZkpXnr65dMqrXb3jhCTa88ETBjj+S/f3gpZZRlfNc+/zBSy0j3lexBGMptnQE2NIRwGU10VzjoLnaQYPbUuqiCSHKXN5BV1GUy4G7VVX9SAHLU1LJtEb3YO5abFcwN4VHRhaLswnFU7zfGeD9zgAOs5Eac5oVlgrq3dZSF00IUYbyDrqqqr6pKIqtkIUptnA8TdfRIBuI0xtJlnUmlpG0+p7+t3vY+NKTGExmMqkkF193I5++/e9GfaxIJILD4chrf+cqZ6HKWG4iyTS+gQht4cO4bSbm1TpZUOek0lG8tTqFEOVtrN3LyZG8yOfzjfEw+QsEAkBuYExvJE1vOIU/ksIfThFNla6reLwEen1cdPWnUa76FOqrv2HA30MkEhn1fuLxeEH3Nx5lLEdH6y0SgSO9Ad7YCZV2E3OqrDRXWXHIKOjTOvo5FaMj9ZafUtZbUa7plirdVjqj8YcDQeKBBAOR47M8GdCZDUzGBsgt9/xk+HHzkgvy3o+mgUHn4oZv/oxkXEcipmflVR8iGddz4H092km97rqTnpgtWcw2DYtVy93bNCy2DCZLtmBlLFcOh+OE5wlgV7/G7v4Y071WFtW7aK5xYDLIOMbjySL2+ZF6y4+kgRwHWhb29sZwOMr3y22w388v/vEv+eJf/5CKypEnzc5HNguxsJ6+LjMBv4lwwEgoYCA8YCQcMJ7wPB4dvxaZxZbB6cng8qZzt8o0FUP3Lm+aiso0VQ1JrI7JdT09S5ZDAzEODcR4fY+euTUOFta7mO6xyqouQkwRkzroTgTHT50pxHVNTYOA34Sv00x/t5m+LjP9PSb6u0z095iJR0YeTM1DrVSzJddSNVuP3QzH7ebky+DZLCQTepIxPYn4iffJuJ5EzEAiZqCv6+xdDQ53muppSaqmJalqSFLdkKR6WoqaGQlszokdkFMZjV3dIXZ1h6iwmlhQ72RRvQu3LFMoxKQ2ltHLy4BmRVGWqqq6o4BlmhTGa3rP0UFK2SwEe430dFjobrPQ026hu8NKT7uFZPzMLXuLLUNlfQpPTSrXsvSkcXrSuLwZnMOP02QI4XI5zriffGkaJKJ6QgEjoX4joYFTbwGfkb5uM5GgkUjQSPsu+yn78dQmaZiVYNqcOA1z4kybnaBqWpICrdZVVIPxFJvbBlDbAkzzWFnc4GJujQOjdD8LMemMZfTydqCxgGURZzULUHjh4VoO7bVxeL+VWPj0rVanN01dU4LqhiSV9Ukq61NU1iepqk9hr8gwkp7M8RrXpNeDzalhcyapnXHmcXiaBoN9Rvq6cq313iNm+rpM9B424ztkIeAzE/CZ2bXJNfwek0WjYVacGfPjzFwYZdbiGN661Ih+33KQJTu8MtKbew3Mr3OyuMFFrUvm/woxWUj38jgZy/SecMBA2y4bna02Du21cWivlWgo96d6/alj77dXpKmfmaBuZoL6oVvdzASOiok/KluvB09NGk9Nmubzoif8LJOBviNmjhyw0nXQypEDFroOWAn2mehotdPRaued5yqB3AnIrEVRZi6KMXNhlBnz4pgs5Tst7KhEOsP2w0G2Hw5S7bSwpMHF/DonVpOMfhZiIpOgW0LhQD8XX/dZ5l/0ed565ml2berjn77ajP/wqS0bhztN4/wYM+bFaZwXY/q8OBWV6QnTiiskgwFqG5PUNiZZfsXg8PbIoIEj+610tNpo32WjfbeN8ICRHe9UsOOditx7jVmaFsSYuzzC3OURmhZEMRbpMmq+g+Z6wwne2Jvgrf19zK1xsmSai+meCT1FXogpS4JukWkadB2wsm+bnWz2aT54y86GF4zAx4ZfY7JoNC2MMXNhjBnzYjTOj+GunpoBdjQcFRnmXRBh3gW5vvFsFnoPm4cDcPsuO91tFg622DnYYufVX9ZgtmrMXhpl7vkR5i2P0DAnPm7Xhcc6aC6jZWntCdHaE8JrN7O4wcWiepesAyzEBCJBd5xls9DTYWHf+3b2f+Bg/weOU67FurwpZi2OMXtJlFmLo0xrjmOQv8yY6XRQMyNJzYwkylVBIDdlav92O/ved7DvfQc9HVZaVSetqhMAR0WaBUqYRSvDLLgoPOpR0t/7zMpzvqYQ+aYHokne3t/HhgP9NNc4WNLgYobXJlOPhChz8tU+DgZ6TOzZ6mDv+7kgGw6cWM3e2iRzz48wZ1mUWUuiVDVMnME+E53NqbH0kjBLLwkDEOwzsm+bYzgIB/wmtv7ew9bfe9Drs8xeGmXRyhCLVoWpmZ4su7+Tls2y1xdmry+M22ZicUMFixtcsv6vEGVKgm4BxCO51tOerU72bnWcck22oipF83kR5p6f68asrJeF7Qst3+ul7qo0F/1xkIv+OEg2C77O3Ijo3ZudHNxxrHfi+f+A6mkJFq0Ms+zSQWYujp22G/p7T206JSMVFCffdDCWYsOBPjYe7Gd2tYOl01w0SutXiGGpjMaB3ijWdOnm+UvQzYOmweG9Vlq3Otmz1UH7Ljta5tgXm9WRyV0jvCA3WKccW0iTTSGSjOh0UNeUpK6pj9Wf6SMW1tO6xcmuTU52q056j1j4wzMW/vBMFS5viqWXhjjv0kFmL4uekCzkdI4Omrv4uht498UnCQ305lXGkdCyWfb7w+z3h6mwmlgyLXft12GRj7uYerLZXCa43T1hDvgjJDMaa+aUbiCifApHKBwwsGerg92qkz1bnEQGj1WdXp9l5uIo8y8MM//CCI3zY+f8EhajM9JkI4Vcn9fm1Fh+xSDLrxhEy0D7bhst77rY/lYF/d1mNjxfyYbnK3FUpFnyoRDLLh2kYW70tPv60nd/PPz407d/d8RlGKvBeIoNB/p59+AAc6rtLG6oYGaltH7F5NcbTrC7O8yenjCRZLrUxRkmQfcMtAx07LEND7I5tNdKNnvsi8pbm2SBkguyc8+PTPi0hOLs9AaYvSTG7CUxPvoVH4f3W9n+Vi4A+w9b2PRbL5t+68XqnM75Hw5xweogs5dGyyZDVjabZb8/wn5/BJfFyKIGF4sbKnBZ5StATB6heJrWnlyg7YskSl2c05JP3HEiQQOtWxzs3uyidYtjOCEFgNGk0bRwPwH/TXz2L/6F2Uvc0mVcRGdroRZ7fV6dDmbMjTNjbpxrvuSnp93CB29V8MFbLnrarWx8ycvGl7y4q1IsXx3kgtWDTGuOl83/SyiRZlPbAJvbAjRV2VjSUMHsKjt6fZkUUIhRiCUz7PXnAm13MHHcanLlaUoHXU2Dw/us7N7sZPdmJ517bCe0ZivrkyxaEWaBEqb5vAjP/8c9HNyxkfdf/zFzlk78Rdcni/G6XjqSwVk6HdTPSlA/y8/VX/RzYKdG66Z63nutggGfmTeeruaNp6upbUwMBeAg1dMKN5BuLKtUZcnS3helvS+Kw2xkYb2TRQ0uvPZJuOalmFSSaY2DfVFae0J09sfQTl51pYxNuaAbDenZsyU3MKZ1i/Ok6TwJ4HXgReAl+rv38vZz8PZzJ+5jpIsTwOiuIYr8jNf10nwGZ9XNjDFnsY9rvuSjfZeN9153s+2NCnydFl55tJZXHq1l9pIoylUBzrtscMzLFxZqlapIMs2WjgBbOgI0uK0saaigucaB2Vgm/eNiyktnNNr7Y+zxhWnrjZI+eVHvCWLSB11Ng0N7rexWc4G2Y7eNrHasNeupTRLwPUwu0P4eOP1AGDH5jOfgLJ0OZi2OMWtxjI+v7Wbv+w7ee83N9rcrhjNiPfNgPUsvCXHRlQHmLY+gP8Pgu5GU82xlHO2JX1cwTlcwzpt79cytdbKkwUW92zqqfQhRCBktS0d/lL2+CAd7cyOPJ7pJG3QPH4bv/p2JZ19YTiRwLLmuwZhl9rIIC5QwC1eEqWtKoNNdAVxx1v0V+7phKeh1OuxmAw6zkYQxTXWlE5NBh8mgx2zQ5+6NOswGPfrjLlAefXj8FcG0liWVyZLKaEO37NA2jWRaI5rMDN8m6hnrSBmMsFCJsFCJ8Kmvd7P9LRfqeg8Htjt473U3773uxl2V4sI/DqJcHTjr6kvFlMxo7OwaZGfXIJV2MwvqXSysc+KUwVdiHGlalkOBXIv2gD9KIj3xF3A53qT99Did8KsnDKTTOtzVKRauCLNQCTN3eQSrffRf8sWcZzkedDodLosRt92Ex2bCZTXiMBuwD90cZiNWk354KonP56O2trYoZUumNSLJ9HAQDifSBGNpgrEUwViKwXia7DhcsynF4CyrXWPF1UFWXB2kv9vElt+5Udd76O8289qT1bz2ZDWzFkdZuWaA8z48iMWWPWM5i3ki2B9NsuFAH+8e6GdGpY2FdU6aaxyYZM1fUQBHW7T7/REO9E6+QHu8SRt03W7493UptgT2MHuhfswjR0s1z3K0bCYD1U4zlQ4zHpsJ99DNZTViKNPRqWajHrPRjPfUteqB3JnvYPxYEB6IpvCHE/SFk+PW3VSMk6zK+hRXfaGXKz/fS1uLjc2vetj2ppu2nXbadtp59t/rOf+KQVZeHaBpYeyU/+FSnAhmydLZH6WzP8rre/TMrXGwoN7FDI9V5v6KUUllNDr6Y+zz567RToau45HQjUcL4nh33XVX9p577hnXY5xJMq3xL7/dcdq0fBOdQa+j0m4eDrDVTjNVDnPBsg4Vs6Wbr2w2F4x7w8nhIOwPJwnFS5dmMxKJjOn/LRHTse1NN5te8dC+89hZSF1TnBVrAlz0x0GcnvJrBbgsRubVOZlX66TWderSlOcyEf7fytFEq7dEKsPBvigHeqN09EdJlSjQrpljY/7MaeO2/7vvvpv77rvvtGehk7alC9Dd3cV//t3X+NO/+dGop1OUEx06PHYTdRWW3M1lpdppLtuWa7HodLrhlnxzzbFAF01m6ArG6R7MDQjyhRJktIkxpcBiy7JyTYCVawL0dJjZ/IqHLb/z0NNh5fmH6nnpP+tY+qFBVl0ToPn8SNkk3wgl0mztCLC1I4DXbmZerYP5dU6ZfiQIxdMc6M0NhDoUiI/LpaKJZFIH3e/fdx8du94f83SKYjMb9DR4rDRUWIeCrAWLSfJKjpTdbKC5xjEciDNaFn8oQddgnK5ggsOBGPFU+bUWT1bXlOT6r/q49mYfOze62PSyh9YtTra96Wbbm24q65OsumYA5aogFZXlk+ZuIJpkU1uSTW0D1DgtzB9qAUv2q6nDF0rQ1hflgD+CP1xemaEG+/18+e/u5Oknf0V9fX3Rj5/3p0BRlL8AfIBbVdUHCleks7NaRz91odzn1VqMBqZ5rExzW5nusVHjNEt2oAIy6HXUu63Uu61c0JjrlvaHk3QOxDg0EONIIF7WI6gNRlh2aYhll4YI+I1sesXD5pe99Hebeem/6nj5kVoWrQpx8bUB5l8YPuPUo1LwhxP4wwne3t9HrctCc42DOdUOKh3SAp5MkmmNzoEYbUPJVsop1/HJ1j/2IFvVzdx7773cf//9RT9+XkFXUZTLgCpVVf9ZUZTvKoqySlXVjQUu26RlMuiZ4bUxw2NjuifXVSyDUIpHp9NR67JQ67JwUZOHjJalKxgfDsI9g+WbSs5Tk+bqL/Ry5Z/0smerk3df8rBro4uWDRW0bKjAU5vMdU9fHcBdXV5ffL5QAl8owYYD/VTazcwZ6o3I5xqwKL2BaJL2vhht/VEOD5R3VqjTzXVft24d69atAyAejxetLPm2dK8Ddg093jn0/IxB1+fz5XmYU3V0dIzodd/5znf45S9/id5oQkunuOiqT3P92m+f832RSGSsRTyFDh1VDiPTKsxMd1uodZqGrscmyMYS+GMFP+SYBQKBUhehqMxAsxOanSZiKQOHAgk6AwkODyZHNdijmB/exsURGhf3EBow8f5r1Wz9XTUDPVZeebSWV39Zw/wLA1x0lZ+5y4Nl1fqF3Oes0z/AG4DTYqDKlGZBIEF9hRmj9PKMWDE/p4m0xpHBJEeCSY4MJgklyuukbiwKGaPOJd+gWw0MDD2OA2ftGC/F6LpQKMQtX70V44LL2fba84QGeos6itluNjKz0kZjpZ0mrw2bucy+9UZgIo2KLLSZ03P3GS3LkUCMg31R2vqiBGPnHhld7NHyDgdcc9MgV39hkH3bHGx80cuODS5aVd2hmf4AABkvSURBVC+tqhdPTYoVawZYeXUAT035fVFmgfZwhN7DSYxdaWZ4bcyssjGr0k6FzXTO90914/U5Pb4HqLM/hi90tAdIB0YLDmNxeijGkl/8qKOXDY/ObTeZzKTTKb761a8WvYs536DrB47OZ3ABfYUpTuE88cQTw1OGmos0r7bKYWF2tZ3ZVXbqKizSZTwJGPQ6GivtNFbauXwe9EeS7PNH2OeLlN3SYXo9zL8gwvwLIoQGDKivetj4Wy99XWZe/UUt6x+rYaESZtW1AyxcES7LNZ/TmkZbX4S2vghvAJV2MzOr7DR6bTS4rZILehylMxrdgwmOBOMcCcToCibKYqxDIfKLZzLQscvG/m1h7BW3Eo/8T/70T/+dnp7uApf23PINui8C1wK/AhYDvy1YiSYQvU7HdK+NOVV2ZlXJWflUUOkws9JhZuUsLwPRJPuHAnC5jdB0eTP80Y19XPGZPvZvc/DuSx5aNlSwa5OLXZtcVFSlWHF1bmpSZV3p5jWfS380SX80yXudAXQ6HXUuCzO8ubEQDW6rZMQag2Rao3swzuFAnMOBWNlMrTtTrvHjB8SOZMBrwG9kz9bceuh73nMQjxiAFwCwWDRuuul+Lr20+L9vXkFXVdW3FUX5I0VRvgwEVFV9s8DlKlsmg56mSjtzaxzMqrLLmfcU5rWbUWaaUWZ6CcZS7PdH2Lo/STldotfrYd4FEeZdECEc6EZd72HTbz34D1v43eM1/P6/q5l3QYRV1wyw+OIQxjI+b8xms3QP5uZfq+25k97aoSB8dGpdoZLDTDaalqUvkqQnlKB7ME7PYIKBaGpSzZnNZKB9l53dm5zs2uyku+3EmS41MxIsuCiXc3/tx7Kct6ChJOXM+z9UVdX/U8iClDOL0cCsKjvNNQ6aKm1ydi1O4baZuLDJwwxrEpPTw15fhD09YQai5bF4AYDTk2H1Z/q44n/0cWC7nY2/9bD9rQr2bHWyZ6sThzvNRR8JsnLNAHVN5VPuM9GOC8JHuSxGao+b317rsky5E+OMlmUgmqQvksIXStAzGMcfSpZFV/FIHN+KPVd+8WhIT6vqZNcmF7tVJ7HwsWsmZqvG3OWRXKBVwlTWH+vRsVptxfllTkNOC8/AYjQwp8bO3BonjV7blM/+JEbOazezclauC9oXSrCnJ8xeX5hwmYz21Omg+bwozedFif7Pbrb8zsPG33roabfy5q+rePPXVcxcFGXlmgDnXx7EYps4raFQIk3IH2a/PwzkZg547SYqnWYq7SYqHWYq7WY8dtOE/0xns1lCiQyR3gh9kSR9kSS94SSBaKqsp++Mxsn5xQcHeulpN7Nzo4tdm5207bSfsFRr9fQEi1aGWbQizOwlUYzm8qsHCbrHMRv0zK52MK/WQVOlfcJ/KEXpHZ0PfGlzJYcDcVp7wuzzhcsmubvdpfHhT/Zz2Sf66Wy1sellD++/UUH7Ljvtu4YWXbg8yMo1AWYuOnXRhXKXJTt8Xfh4Op0Oj81EpcOE127GaTHishpxDd2XS+s4ldEYjKUJxlPD98FY7vFgPMVgKIzDMXnXAP/Sd39MJg0HW+wYTVfQfdDFD//sWGIVvSFL8/kRFq0KsWhFmJoyWRbzbKZ80DXq9cyssjO/1sHMKrt0HYtxodPpcglRvDYun1fFgd4ou7pDHOqPlUUiDp0OmhbGaFoY4+Nf62bbHyrY/LKXgy12Nr/iZfMrXmpmJLjoI0Eu+kh5Tj0ajWw21wWb6/4/dW6+xWjAYTEMLYFpxGLU524mPRajAYtRj9Wox2zUY9Tr0Ot06HS5v7Nex7HngJbNBc+0liWdyZLWjq0vnc5oRFMZYskMsZRGdGiJy6PPJ/MSd2cTC+e6jVs2umg9qdvYXpFm0Yowi1eFmHdhBJujPE5gR2pKBl2dTkej18aCOidzqh1lc1YrpgaTQc+COicL6pyE42l294Rp7Q6d0horFbM1y4qrgqy4KojvkJnNL3tQf+fBf8jCb39ey8uP1NB8fgTlyiDLLh3EbC39SUOhJdIZEukM/ZHy+JtMBf3dJlredbFzo5MD2x1omWPdKrWNCRavCrF4VYiZi2Jll+xlNKZU0K11WVhQ52JerUNGOYqy4LQaUWZ6UGZ66BmMs7MrxF5fpGQtnJMTEdTOSPLRW3xcc7OPPVucqOvdtGxwse99J/ved/Kbn9Sz7LIQypUBZi+Nls2qR6L8aRp0ttrYudHJzo2uE0Yb6/VZ5iyLsPjiEItXhamZPnlOfiZ95HFaDFw008uCOqckWRdlra7CSl2FlcvmVrHfH2FnV4gjgXhRu5/PlIjAYCA3QGVlmGhIz7Y33ajr3XTstqO+6kF91YO7OsXy1UEuuGKQac3x4eu/hcgoJCaHZFzHnq1Odm1ysnOTi/DAsRBktWdYoOS6jReuCGN3Taxu45Ga1EHXbNRzw3nV1NVVlrooQoyYyaBnYb2LhfUuAtEUu7pD7O4OFXT085kSEBx1tpW5fvBSC5d8dIBLPjqA75CZLevdbH3NTcBn5o2nqnnjqWpqGxO5ALw6yJu/HntGITFxDfSY2LU515rdv81OOnWsO8Rbm8y1Zi8OM2dppKzniRfKpA66gKRiFBOax27ikjmVrJrlpWMgRsuRQQ72RcsmqUHtjCTX3uxnzZ/6+fZHbwU+D9yIr7OGVx6dySuPHptDe6ZAXorlNMX40TLQ0WrLtWZP6jbW6bI0LYyyeGWYRReHaJiVmHAj4sdq0gddISYDvV7HrKF0o+FEmt3dIVqOhBiM55fC8XSB7lyJCM5ePoB3hm53AB8Bvg+8RS71XpRcuvbLgWbgZWBfXmWfCKZal3o4YGC3mku52LrlxNHGFluGeRdGWDw0rcfpmZojso+SoCvEBOO0GFFmermoycOhgRg7ukIc8EfGnBDh5EQEoYHeUb3/dIH8yX99h02vxNBhJZuNA7OBBwConxVn2aUhug4O4qoZ+ZKaEyGgFSJJfznTMtCxx0brZie7VSeH9h7N8NQFfBxv7aMsXuVk0aoQzcvKM0lFqUjQFWKC0umOrYAUS2bY1R2i5cgggREsP3g6X/ruj4cff7pAK3NFQ31cMhTI33nuSboOHqR6emC427G7zcqrv6zBWz+dJasiLFwRZs6yKKazfEmXY0AbSZL+oyZid3o2C71HzOx9z8He9xzs/8BxQmvWaNKYc16UZOyvadv1BxauuJNP3lYef5tyI0FXiEnAZjZwYZOHC4davy1dg+z3R0q+aszxgfwz3zwayI+QTsG+bQ62v11ByzsuBrqtvPWslbeercJk0Zi3PBeAf/3AKqDztPuWa8TjKxwwsG+bYzjQDvhOnP1R1ZCkr2sd8CLp1Evs2XLsZyf/beRvcowEXSEmmaOZr2LJDLt7ctd+y2nhBQCjCRYqERYqEf7H7V20vq+jfUc1uzY5OXLAxs6NLnZudAEdwHbgFeAl4A2g/LJhjSZJf7kK9hk5uN3OgR12DrbYT1mlx16RZt7ySO52QYTK+hR/de03SlTaiUuCrhCTlM1s4IJGDxc0ejgciNFyJMQ+f7jkrd+T6Q3QtDDMoouyXPMlP8FeI61bcnM5977nIBFbBiwD/gL4M+Ah9HoTWjbJiqs/y413FKYrvFDGem28GLJZ6OsycbDFPhRoHfR1ndiSNZo1Zi+J5oLshRGmzYmfkvxkMpxsFJsEXSGmgOkeG9M9Ni5PVbG7O0xL12DZpjh0V6dZuSbAyjUB0iloa7Gzb5uDfR84aN/pA/4MTVsLrGPzy130tM1i1pIYjfNjNC2I4a1LlXQaynhcGx+r0ICBzj02OlttdO610dlqJRo68evfYsswa0mMOUsjzF4apXFefFQDoCbCyUY5kKArxBRiNRlY3uhmeaObrmCcliOD7PVFynatVaMJ5i6PMnd5FPCTjP8DB1vs7P/Awf4PfsChPTY6WnV0tNqH3+OoSNM4P0bjgnjufn5sykxTyWYh4DPR3W6hq83C4b02OvZYCfhOzcbn9KSZtTjKnKVRZi+N0jAnjmEMOY3L8WTjZEeXejQZSndWJkFXiCmqwW2lwW3lw3MztPaEaekK0RtOlLpYZ2W2ZllwUYQFF+WmGMWjetp2DrXg9tjoaLURCRrZrbrYrbqG3+f0pqlvSlDXlKB26L6uKTFhg7GmwWCfkc59LgI9HrrbLXS3Wehpt5CInRo5LbYMM+bFmTE/RtP8GI0LYnhq0pM6MYUOHW6bkdoKy/ASmzVOC2ajHp/PV7JySdAVYoqzmAycN8PNeTPc9Awm2Nk1WNJFF0bDateGB2RBrqU34DPR2ZoLwIf2WDm0z0Z4wMi+ASP7tjlOeL+jIk1tYxJPbRJvbRp3dQpvbQpPTQpPbapky8alUxAJGhnsN9LfY6K/20x3Wx+7Nq3F5vwlwd5GMunTry7h9KapnxmnfmaCaXPiNC6IUzsjMaFX5hkJk0FPvdvKtKGTyVqnGYup/H5pCbpCiGF1FRbqKmqGF11o6QpxJBArdbFGTKeDyroUlXUpzr98EMi1CoN+Ez0dFrrbLfg6zPR0WOjpsBAZNHKwxQgt9tPuz2LLUFGZxubMYHNpuXtnBvvQvc2pYbZo6AxZ9HowGLLoDbnF1fVD29JpHemEjlRST+rofVJHKqEnGdcRGTQSDhgIB41EggbCASPxyOmCxT3ABmLh7wM/xelN462LMX1OivqZCepnTezW+2hZTQamua1M81iZ5rZR4zSj15d/0z2voKsoig34FqCpqnpfYYskhCi14xddGIgm2dUVYnd3mEiy/KbrnIteD966FN66FAtXhIe3Z7MQ7DXiP2Qh0Gsk4Dflbj7T8ONEzID/cClaSxnAD/QAK4DjE548CDxIeADCA9C5e2rMgzUb9Ez32mgcunntpgmZWz+voKuqakxRFBX4UIHLI4QoM167mQ81V3Hx7Eo6BmLs7g5xoLf0iTfGSqcDT00aT83pTySyWYiF9YQGjMTCBmJhA9GwgVhIf8LzdFKHltGhaQzd69AyQ48zOgzGLCaLhsmcxTh0bzJrmCy57Y6KDE5PGqc7g8OdxunJ8L3PLoIiLOlYzik19Tod9W7rcJCtc1kmREv2XMbSvTzi+QalvGgdCARKduyJTOotP5O93uzAhTU6lnjtHOyPs7c3hj+cX9rJ48Xj8XO/qBT04KzK3Yrpe09tPOH5c+vuY+urv8FgNJFJp7joqk9z/dpvE4/HsVqtRCIjz119vJd+/m+0tWzhpZ//G9ev/XYhij4mFVYjM9xmprst1LtMmAx6IAWJFL0FHONXys/pOYOuoijfAeaftPkZYMSlrq2tHWWxCqvUx5+opN7yM1XqrXFabs2g/kiS1p4wrT1hQnmuegTgcDjO/aIpKhEePGUO7NH6Gmm9nW0NZfWVp1FfeXr4ebG6q00GPTO8NpoqbcystOO2FW9B3VJ9Ts8ZdFVVvfd02xVFWV3w0gghJpxKh5lL5lRy8WwvXcEEe3xh9vrCxFNTY0BPMUyEObAjVekwM6vKTlOlnWluK4ZJ0GU8GjJ6WQhREDqdLjeS1GPlw3Or6OiPsscX5mBvlFSmPJNvTCUnt16LlbbRbNDTWGmjqdLOzEo7LuvUDjv5jl42khtEtURRFK+qqgOFLZYQYiIz6HXMrnYwu9pBMq3R1hdlnz9Ce1+0bLNfTTXjmbZxqrdmzybf0ctp4LTdzkIIcTyzUc/8Oifz65ykMhod/TH2+cO09UZJSgu4ZArZZW0xGmistDGz0kaT145zirdmz0ZqRghRNCaDnuYaB801DtIZjc6BGPv9Eba3lzYBRzlPnSlHOp2OOpdlKNDaJ810nmKQoCuEKAmjQT/cBb2sMotmddPWF6WtL1r0HNDrH3uQtpYtrH/sQVmO7gy8djONQ2s1z/BYyzLF4kQgQVcIUXI6nW54AYZL5lQSjqdp64/S3helcyBW0IFYZ5s6s+GFJ9jwwhMnbJsK2Z5Ox2E2MuO4DFDSZVwYUotCiLLjtBpZOq2CpdMqSGc0ugcTHA7EOBSI0x2Mo2UndjascuS2mZjmsTLdbaPBbcVjL96c2alEgq4QoqwZhxIozPDaWAWkMhrdwTiHAnEODcTwhRKjCsKlmjpTTnQ6HdUOMw3DCwZYcVgkHBSD1LIQYkIxGfQ0VtpprMytDJRMa/jDCXoGE/SEcvejyYw1nlNnyoEOHR67ibqhdWXrXBaqnWaMhtMvDSjGlwRdIcSEZjbqme6xMd1jG94WTWboGYzTE0rgCyXoDycJJzJkT7OIwGTK9mQ26Kl0mPE6TFQ5zFQ7c0HWbJQAWy4k6AohJh272TA8MvqoZFpjIJqkP5Ki/+h9JEkonj5tMC5XOnTYzQZcViM6m8bsaVVUOcxU2k25bRNwubupRIKuEGJKMBv11FVYqauwnrA9o2UJxdOEE7lbKJ4mNHQfjqeJJDMk01pRArNep8NqMmAzGbCZ9dhMBiqsRlxWE26bEZfViMtiHO4a9vl81NZ6xr1conAk6AohpjSDPnfN82yjdbPZLMm0RjytkUhrJFIZEkPP05ksWjaLls297ujjo/cGXe4YBr0eo16H0aDDqNdh0Ofuh4OsSS9zX6cACbpCCHEOOp0Oi8kgQVGMmVxdF0IIIYpEgq4QQghRJBJ0hRBCiCKRoCuEEEIUiQRdIYQQokgk6AohhBBFIkFXCCGEKBIJukIIIUSRjDo5hqIoHwG+B0wHvq6q6kuFLpQQQggxGeXT0q1QVfXDwK3ADwtcHiGEEGLSGnXQVVX1N0MPNwNdhS2OEEIIMXmdtXtZUZTvAPNP2vyMqqrPANcB943kID6fL7/SFUAgECjZsScyqbf8SL3lR+otP1Jv+SllvZ016Kqqeu/ptiuKUg04VFX975EcpLa2No+iFU6pjz9RSb3lR+otP1Jv+ZF6y0+p6m3U3cuKojiA61RVfVhRFKOiKFXjUC4hhBBi0hnV6GVFUSzAC4BLUZQ/B9zAheNRMCGEEGKyGVXQVVU1Aawen6IIIYQQk5skxxBCCCGKRIKuEEIIUSQSdIUQQogikaArhBBCFIkEXSGEEKJIJOgKIYQQRSJBVwghhCgSCbpCCCFEkUjQFUIIIYpEgq4QQghRJBJ0hRBCiCKRoCuEEEIUiQRdIYQQokgk6AohhBBFIkFXCCGEKBIJukIIIUSRSNAVQgghikSCrhBCCFEkEnSFEEKIIjHm8yZFUdYAtwP1wPWqqvYUtFRCCCHEJJRvS7dDVdWPAb8BVhWwPEIIIcSklVfQVVV119DDELC+cMURQgghJq9zdi8rivIdYP5Jm58Zeu93gM6h52d0991351s+IYQQYtLQZbPZvN+sKMpK4G5VVT9auCIJIYQQk9NYRy8fAFoKURAhhBBishv16GVFUXTAi8DvgCDwvwtdKCGEEGIyGlP3shBCCCFGTpJjCCGEEEUiQVcIIYQokrwyUpUrRVH+AvABblVVHzhu+3zgs0AUeE5V1T0lKmJZOku9fQ64A6gAblJVVS1REcvSmertuJ8/DDyqqurrxS5bOTtbvSmKshD4MLBDVdUNpShfuTrL5/RTQNXQ06iqqo+VonzlSlGUy8nNsvnISds/BFxKrvH5n6qq+opRnknT0lUU5TKgSlXVRwGvoijHZ8r6MfAj4AHgvlKUr1ydqd6GBsxFVVVdBfwQuKeExSw75/h/Q1GUjwHOkhSujJ2t3hRFWQDcqqrqQxJwT3SO/7dvqqr6H6qq/gfwldKUsHypqvomYDvNj/6R3Hfb4xTx+23SBF3gOuBopqydQ89RFMUGNKuqGlZVNQHMVhRlUrXwx+i09aaqalZV1WeHtm8GukpQtnJ22noDUBRlNrlepF2ned9Ud8Z6A+4H2hVF+fFQkBHHnK3etiiK8veKoijAT4tesokhefyTod7P9ND3XAe53pWimExBtxoYGHocJ7cYA4AXGDzudWmgpojlKndnqrfjXQn8S9FKNDGctt6GTuiuVVX1N6UqWJk7U705gFnkeqP+GXhSURRzKQpYps72Of0u0Az8AHizyOWaqI6vT8jFiaKYTEHXD9iHHruAvqHHfYD1uNfZgUARy1XuzlRvACiKMhdoV1V1Z7ELVubOVG+XA19UFOV14GbgXxVFmV700pWvM9WbGYipqqoNtTyOcPoTwKnqbJ/TfwD+jFxX6X8XuVwT1fH1CZAo1oEnU9B9EThv6PFi4GVFUdxDXcrtiqLYFUWxAp2qqsZKVsryc9p6A1AUpQ44X1XVpxVFcQ61RkTOmf7ffq+q6odUVV0N/Bdwh6qqh0tUxnJ0pnobABKKohy9Du4HpN6OOePnFDhPVdWQqqovAKaSlG6CUBTFoCiKS1XVvQw1xhRFmQO8XqwyTJqgq6rq20BcUZQvk2vJBoB/H/rxt4E7gW8B/6s0JSxPZ6o3RVGqgJeBv1YURQXeIDf6W3DO/zdxBueot9uBu4dGzX9fVdVMiYpZds5Rbz9SFOUbiqJ8GvhZqcpYrhRFWQY0K4qylNy18L8d+tH/URTlTuCLwN8UqzySkUoIIYQokknT0hVCCCHKnQRdIYQQokgk6AohhBBFIkFXCCGEKBIJukIIIUSRSNAVQgghikSCrhBCCFEkEnSFEEKIIvn/pgHNs2Ghx1oAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 576x216 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Get into evaluation (predictive posterior) mode\n",
    "model.eval()\n",
    "likelihood.eval()\n",
    "\n",
    "# Test points are regularly spaced along [0,1]\n",
    "# Make predictions by feeding model through likelihood\n",
    "with torch.no_grad(), gpytorch.settings.fast_pred_var():\n",
    "    test_x = torch.linspace(0, 1, 51)\n",
    "    observed_pred = likelihood(model(test_x))\n",
    "\n",
    "with torch.no_grad():\n",
    "    # Initialize plot\n",
    "    f, ax = plt.subplots(1, 1, figsize=(8, 3))\n",
    "\n",
    "    # Get upper and lower confidence bounds\n",
    "    lower, upper = observed_pred.confidence_region()\n",
    "    # Plot training data as black stars\n",
    "    ax.errorbar(train_x_mean.numpy(), train_y.numpy(), xerr=train_x_stdv, fmt='k*')\n",
    "    # Plot predictive means as blue line\n",
    "    ax.plot(test_x.numpy(), observed_pred.mean.numpy(), 'b')\n",
    "    # Shade between the lower and upper confidence bounds\n",
    "    ax.fill_between(test_x.numpy(), lower.numpy(), upper.numpy(), alpha=0.5)\n",
    "    ax.set_ylim([-3, 3])\n",
    "    ax.legend(['Observed Data', 'Mean', 'Confidence'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is a toy example, but it can be useful in practice for more complex datasets where features are more likely to be missing."
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
